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I. INTRODUCTION 



We investigate a model of globally coupled conservative oscillators. Two different algebraic poten- 
tials are considered that display in the canonical ensemble either a second (0'*) or both a second and 
^ ' a first order phase transition separated by tricritical points {(f>^). The stability of highly clustered 

^T^' states appearing in the low temperature/energy region is studied both analytically and numerically 

, for the ^''-model. Moreover, long-lived out-of-equilibrium states appear close to the second order 

QQ ■ phase transition when starting with "water-bag" initial conditions, in analogy with what has been 

found for the Hamiltonian Mean Field (HMF) model. The microcanonical simulations of the (f)^- 
model show strong hysteretic effects and metastability near the first-order phase transition and a 
• narrow region of negative specific heat. 

o 

B 

I ; I The treatment of long-range interacting systems remains a challenging issue in thermodynamics and statistical 
. mechanics Serious theoretical difficulties arise because internal energy, entropy and other thermodynamic quan- 
tities are no longer additive, i.e. a part of a system has not the same thermodynamic properties of the whole. This 
I , originates unusual effects, like negative specific heat and the inequivalence of statistical ensembles even in the limit 
' ^ ' of infinite number of particles (See Ref. .23 for a recent revievsr emphasizing different examples such as gravitation, 
plasmas, fluid mechanics,. . . ). Relevant physical examples displaying such anomalies are known in Newtonian gravity 
but also in plasma physics (although in the latter case the screening of attractive and repulsive Coulomb interactions 
may mitigate them). 

As usual in theoretical physics, the study of simple toy models proves to be of major importance to attack more 
complex and realistic systems. In particular, simple mean-field models with infinite-range interactions turned out to 
be extremely useful. In spite of the fact that they are constantly used in statistical mechanics to describe cooperative 
' phenomena, it is somehow singular that violation of additivity has hardly been recognized in the past. A reason for that 
T-H , lies perhaps in the fact that the thermodynamic limit is performed resorting to saddle-point techniques, which puts 
' the Hamiltonian in the explicitly decoupled form, thus hiding the difficulties inherent in the long-range interaction. 

Indeed, ensemble inequivalence (for example between microcanonical and canonical ensemble) has been observed, 
. producing effects like negative specific heats, which are the counterparts of the ones known in the gravitational 
context 1]. 

The advantage of such models is that their canonical thermodynamics can be exactly derived by performing the 
mean- field limit (the infinite iV limit at fixed volume), which is a reasonable surrogate of the thermodynamic limit 
(the infinite N limit at fixed density). Contrary to the usual belief, an exact microcanonical solution is also feasible 
'"^ • for such non trivial Hamiltonians, using large deviations techniques but the results will be presented elsewhere 0. 
' Here, for what concerns the microcanonical ensemble, we will mainly limit ourselves to show the result of numerical 
O • simulations, which, because of the mean-field nature of the interaction, require only 0{N) codes (instead of the usual 
y, • 0{N'^)). Moreover, further insight can be gained from solving the one-dimensional collisionless Boltzmann-Poisson 
^ , equation for the single-particle distribution function, which becomes exact in the N ^ oo limit (at all finite times) j^. 

In the present paper, we investigate, both analytically and numerically, two simple mean-field models which we 
denote as "0'*" and "0^", which display respectively second (0'*) and first and second order phase transitions separated 
by tricritical points {(f)^) in the canonical ensemble. Of the former model, we investigate in addition the dynamical 
formation of clustered states at low temperatures and we study their destabilization. The presence of quasi-stationary 
out-of-equilibrium states is moreover revealed close to the second order phase transition, in analogy with what is 
found for the Hamiltonian-Mean-Field (HMF) model (see for a recent review). Concerning the 0^-model, we study 
the phase diagram in the canonical ensemble and we report numerical simulations of hysteretic effects near first-order 
phase transitions. We point out the existence of a narrow region of negative specific heat. 
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II. THE MEAN-FIELD MODEL 



Let us first consider tlie following Hamiltonian 

N 
i=l 

where pi is the conjugate momentum of the variable qi, which defines the position of the i-th particle on a line. This is 
a mean field model since all particles are connected to all others, and the summation in the last term is not restricted 
to neighboring particles. Notice that positive (resp. negative) values of the parameter 9 correspond to attractive 
(resp. repulsive) mean-field interactions. All variables are adimensional and, for the sake of comparison, we have used 
the same parametrization introduced in Ref. (which can be shown to be minimal by conveniently rescaling the 
variables and time). The local potential displays a double well for 9 < 1 and a single well otherwise. The ground-state 
energy per particle is eo = — 1/4 for positive 9 (all particles in a single cluster) and cq = —1/4 + 9/2 in the repulsive 
case (double cluster). 
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A. Dynamics of the Magnetization: The generation of clusters. 

Introducing the time-dependent magnetization 

1 ^ 

^^=]^E* ' (2) 

1=1 

we are therefore interested in the following equations of motion 

q^^{\-9)q,-(^l + 9M . (3) 

We study the dynamics of particle released with a water bag initial condition where positions and momenta are 
uniformly distributed at random in the intervals [50 — 'Wq/'^,1o + Wq/'2] and [—Wp/2,-\''Wp/2]^ respectively. We have 
adopted the symplectic 6th-order Yoshida's algorithm 0, with a time step dt = 0.05, which allows us to obtain an 
energy conservation with a relative accuracy /^E/E ranging from 10~^ to 10"^^°. 

Fig. n shows the result: a coherent oscillating cluster self-consistently moving in the self-generated potential. The 
data are obtained for an initial condition with a small velocity dispersion, i.e. qo = 1-1, Wq = 0.05, Wp = 0.0001. 
Besides the oscillation of the center, the particles display a rotating motion around it, which creates a spiral structure 
(see Figs. 121), as frequently encountered in long range systems; we have found this coherent behavior for a very large 
collection of initial states. Notice that the spiral structure in the center is responsible for the very large peaks in 
the single particle density (right panels in Figs. A similar phenomenon has been described successfully for the 
antiferromagnetic HMF model in terms of shock waves [T(l flll | by considering the associated Vlasov equation valid 
in the N 00 limit. 

We will therefore rely on a Vlasov-like approach. Denoting by f{q,p,t) the one particle distribution function, we 
have here 

T^+PTT+t^ (1-^) / / dafia,u,t)a ^0 . (4) 

dt dq dp [ J 

Introducing a density field p and a velocity field v, as follows 

r + oo 

p{q,t) = / f{q.p,t)dp (5) 



/+00 
pf{q,p,t)dp (6) 
-00 

and neglecting velocity dispersion, we have recently shown |lOj | how to reduce this problem to appropriate hydrody- 
namical equations. A short-time analysis, performed for the HMF model, led us finally to a dissipativeless spatially 
forced Burgers equation. We expect that a similar treatment can be developed for the current model and that similar 
techniques could be applied. A well known property of the Burgers equation without viscosity, is that the solution 
becomes multi-stream after a finite time: the appearance of shock waves in the velocity profile corresponds indeed to 
singular points in the density profile (see Figs.[2Jl. In the original discrete model, this phenomenon would correspond 
to particle crossing; after some time, fast particles will eventually catch slow ones downstream creating the so-called 
spiral dynamics exemplified in the left panels of Figs. 



FIG. 1: Dynamics of the cluster. Evolution of the density p{q) of formula Q in grey scale for short times. The darker the 
grey, the bigger the density. Space is horizontal, whereas the vertical downward direction corresponds to time. One notices the 
periodic motion with the characteristic time scale ujj,^, defined in the text. In this simulation A*' — 4096, 6 = 0.5. 
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FIG. 2: Phase space snapshots of the cluster and the corresponding density p{q) of formula Q at two different times. In this 
simulation iV = 1024, e = -0.2387,6' = 0.5 



B. Stability analysis 

To understand the origin of this cluster and its stability , we will first consider the simplest case of a fully clustered 
state where qi = M. In this simple case, the collective motion is ruled by the equation 

M = M - , (7) 

which can be easily solved using elliptic functions 12]. Integrating Eq. 0, between the initial time, when the cluster 
is released without kinetic energy at the position q = a > 1, and time t, we get 

M = adn(^-^,fc^ , (8) 

where dn is the elliptic delta amplitude function and k = \/2 — the modulus of Jacobi elliptic functions. We 
remind that this solution is periodic, with the amplitude-dependent period given in terms of the complete elliptic 
integral of the first kind 2K{k); the magnetization M will thus oscillate with a frequency ljm — na/^/2K, which will 
be the main timescale of the problem. One notices immediately that the modulus k and the frequency lum are both 
functions of the same parameter, namely the amplitude a, related to the energy per particle e = E/N through the 
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FIG. 3: Critical energy as a function of the parameter Q. The sohd hne corresponds to the resuhs given by the stabihty 
cliarts derived analyticaUy (or, ahernatively, using the Floquet analysis), whereas the diamonds correspond to the results of 
microcanonical simulations for an ensemble of A'^ = 1024 particles. 



relation = 1 + VI + 4e. This solution is interesting in its own, since explicit analytical solutions are not common for 
nonlinear non-integrable systems of oscillators, but one should of course study its stability in order to understand why 
this coherent oscillating cluster emerges spontaneously. Using the equations of motion Q for the qi and introducing 
= Qi ^ M , we obtain up to first order 



Introducing the new variable u = at/^/2, we obtain the Lame equation in its canonical form 



+ [a - v{v + l)fc^ sn^ (u, k)] £,i 







(9) 



(10) 



with a = 6 + 2(0 — l)/a^ and i/ = 2. For integer values of v, many rigorous results are known 0, and in particular 
it is established that there are only i/+ 1 instability regions in the (a, k) plane. The stability charts could be explicitly 
constructed by observing that Eq. itTUI) has the following five periodic solutions 



y = 1 — jsn^{u, k) 

y = cn(M, fc)dn(u, k) 

y = sn(M, fc)dn(u, fc) 

y = sn(u, fc)cn(u, fc) 



with 


a 


= 2 


1 + fc 


with 


a 


= 1 




with 


a 


= 1 


f 4fc2 


with 


a 


= 4 





k^ + l) 



1/2 



(11) 

(12) 

(13) 
(14) 



Thus the above curves a — a{k'^) define the boundary curves of the three {v + 1) non-degenerate instability regions. 
Theses curves are presented in the plane {9, e) in Fig. 

One can also investigate the linear stability of this cluster solution with a standard Floquet analysis, i.e. computing 
the eigenvalues of the 2N x 2N matrix of the tangent map. Here, contrary to usual lattice systems with coupling 
between neighbors, Eq. shows that we obtain N identical second order equations; this is a direct consequence 
of the mean field character of Hamiltonian Consequently, we obtain two different N times degenerate Floquet 
eigenvalues and the periodic solution is linearly stable when the eigenvalues lie on the unit circle in the complex plane. 

At this stage, one derives numerically the linear stability threshold by considering the numerical evolution of two 
different initial conditions (1,0) and (0,1) for the vector (i^, ^). The dynamics is solved by a standard 4th order Runge- 
Kutta algorithm for the time integration of Eq. @, where the magnetization M is either directly integrated using 
Eq. lO or implemented with the help of Eq. ((SJ. For a given value of 9, an energy threshold exists, above which the 
largest Floquet multiplier is greater than unity, and therefore the solution is unstable. The solid line in Fig. O shows 
the evolution of this threshold as a function of the parameter 9. The analytical calculations were directly compared 
with the numerical thresholds obtained by considering a water bag with very small but finite width, i.e. Wq <ti qo and 
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uip <C 1, to make a direct comparison with the above analytical results. Checking on Fig. |31 one gets, apart from a 
slight underestimate, a good agreement between numerics and theory. It should however be remarked that, for the 
finite N systems, stability persists only for a finite time, which presumably diverges as N increases, as it happens for 
the HMF model Ullll. 



C. Equilibrium Statistical Mechanics 



The partition function can be computed by means of a standard Hubbard-Stratonovich transformation. Indeed, for 
a Hamiltonian of the general form 

2 



N 



Pi 



2N 



' N 



the partition function is 



Z = 



/ n dp, dqi e-f'" = ZkZv = {2T:/I3f'^ Z^ 



where the configurational partition function is 
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We use at this point the Hubbard-Stratonovich trick, i.e. we consider the identity 

7^ 



dy e 



Defining 



after some algebra one gets 



ip{x,l3) = In 



+ 00 



dq e 
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(20) 
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where is the configurational Landau free energy. In the thermodynamic limit, one can evaluate the above integral 
by means of the saddle point approximation. The saddle point is determined by the condition x — f39 ^px {x, (3) , (where 
ipx denotes the derivative with respect to x) and can be evaluated numerically in a self-consistent manner. Notice 
that X = is always a solution if the potential V is even. 

Finally, we can thus express the configurational partition function as 

7^2 



Zv 



1 



y/Wxxix, 13) -1 



exp 



A^(^(^,/3)-^ 



(22) 



Up to terms of order O (l/N), the relevant equilibrium observables can be expressed accordingly as a function of x, 
using the following formulae: 
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In the disordered phase, when M = 0, the system reduces to an ensemble of independent anharmonic oscillators. 

In the ferromagnetic case {9 > 0), as presented in Fig.0] the model displays a second-order transition in the canonical 
ensemble, in full agreement with previous results based on the Fokker-Planck approach • The magnetization vanishes 
as (Tg — T) 2 in the subcritical regime and the specific heat has a finite jump at T^. Conversely, in the antiferromagnetic 
case {9 < 0), no transition occurs, and a; = is the only solution of the consistency equation for any value of the 
temperature. The free energy and the internal energy are always given by the above formulae, with zero magnetization. 

As this behavior is clearly reminiscent of the HMF model, that we have already studied in the past 0, 
[T^ . and where the presence of long-lived out-of-equilibrium states was surprisingly discovered, it was natural to 
suspect that they also appear in the present model. We have therefore performed two types of numerical simulations: 
microcanonical ones with a symplectic algorithm (sixth order Yoshida or fourth-order McLachlan-Atela il8]) and 
canonical ones with a Nose- Hoover thermostat 19] (fourth-order Runge-Kutta algorithm). No appreciable deviations 
are observed between the two types of simulations for initial conditions close to equilibrium, see Fig. ^ 
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FIG. 4: Comparison of ensembles for the model, 6 = 0.5, = 512: caloric (panel a) curve and magnetization (panel 
b). Squares and triangles refer to microcanonical and canonical simulations, respectively, while the solid lines are the exact 
canonical solutions given by Eqs. 12511 and I24II . The critical point is located at Tc — 0.264 {ec = 0.132). In both cases the 
initial conditions were qi{0) = 1 and Pi{0) chosen randomly with a gaussian distribution. 

However, a region with clear differences is found for "water-bag" initial conditions: see an example for go = in 
Fig. El This is strongly reminiscent of similar observations made on the HMF model 0, |23| . The fact that some 
points lie on the branch with vanishing magnetization also in the subcritical region (see Fig. [SJs) indicates that this is 
a metastable state in the microcanonical ensemble. On the contrary, let us notice that triangles below the theoretical 
curve in Fig. are due to finite size effects and would disappear for larger N values. A careful study of the numerical 
results for very large integration times shows a systematic tendency of these points to converge towards the equilibrium 
state indicated in Fig. 0] by the solid line. This attests the metastable character of these states. 

Series of microcanonical runs for the repulsive case have shown that metastable states may possibly exist also 
in this case, and we suspect that they may be related to the existence of a stable cluster in the energy region 
— 1/4 < e —0.097. This metastability is thus of dynamical rather than of thermodynamical origin. 

Summarizing, the 0* model has emphasized the striking appearence of a cluster and also dynamical differences 
between microcanonical and canonical ensembles, presumably related to slow relaxation towards the final Boltzmann- 
Gibbs equilibrium state; ensemble inequivalence appears only in a transient. Indeed, it has been recently reported 
in spin systems gH, that true ensemble inequivalence occurs in regions of first order phase transitions. It would be 
therefore very interesting to exhibit a dynamical mean field model of the polynomial class with a first order phase 
transition. This is the purpose of the next section. 



III. THE MEAN-FIELD 4>^ MODEL 



The simplest generalization of the previous model is 
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where D and r are two independent parameters (also in this case it can be shown that this parametrization is minimal). 
The main interest of model H2t)|l lies in the fact that it may exhibit a first-order phase transition for a proper choice 
of the parameters, and therefore possibly ensemble inequivalence. This can be realized by first considering the zero 
temperature limit, where equilibrium states are given by the minima of the function 

r — D ^ ,„„, 

^eff-^-'-T + y • (2^) 

For < r — D < 1/4, such polynomial admits three minima located at a; = and x = ±x+ and two maxima at 
X — ±X- where 



xl = -±-Vl-^r-D) . (28) 

A first-order transition can thus be expected within this parameter region. Furthermore, in order to have a first-order 
phase transition at T = we must impose that the two minima attain the same value (equal to zero). This conditions 
holds for r — D — 3/16 and we can at least hope that close to this parameter values the transition persists also at 
nonzero temperature. We checked that this is indeed the case by computing the free energy in a self-consistent way, 
as explained in the previous section. The transition exists in a very narrow region below r — D = 3/16. 

It is useful to consider the expansion up to sixth order of the Landau free energy of the (j)^-model in order to 
determine the critical line of second order transitions and the tricritical point. We obtain 

'^+0{x'), (29) 



(30) 
(31) 



(32) 

The numerical solution of a = yields the critical line of second order transitions, see Fig. The tricritical point, 
separating first and second order phase transition, is determined by the more restrictive condition a = b = 0. Tabic ^ 
presents some values of the tricritical point as a function of the parameter D. 



X dx bx 

0fL {x,l3) = -ip{x,f3) = const. + — + — ■ 

where 

a{f3,r,D) = ^-{q') , 
6 

with 

Jq^exp{~pViq))dq 



[q ) = 



Jexp{-pV{q))dq 




FIG. 6: Panel (a) shows the phase diagram of the <f)^ model for different values of the coupling constant D. The solid 
(respectively dashed) line marks the second (resp. first) order critical line, and the full dots the tricritical points. Panels (b) 
present the Landau free energy at a first order transition close to the tricritical point Ttr and at a second order transition 
occurring on the line a = Q, Tc = 0.205. 



The canonical thermodynamics in the case of a first-order transition is further illustrated in Fig.[7| Three branches 
of solutions (two stable and one unstable) exist from T = up to T = T' where a saddle-node bifurcation occurs. 
The m = branch is stable at all temperatures, while the stable (upper branch in Fig.0) and unstable (lower branch 
in Fig.[71l TO 7^ solutions meet and collide at T = T' . Notice that this is at variance with the Blume- Emery-Griffiths 
model 21\ (BEG), where the three branches do not extend down to zero temperature. 

We have performed some simulations in the canonical ensemble to check this caloric curve. The results are in 
agreement with the theory and, as expected, display a marked metastability around the transition point (hysteretic 
effects). More specifically, three different initial conditions were adopted: (i) all qi = (ii) all qi — x+ (iii) random 
distribution between qi — and qi — x+. In all cases, the pi were initially chosen according to a random gaussian 
distribution. Some microcanonical data are reported in Fig. [T] for an initial condition of type (ii). 

A peculiarity of this model appears in some region of the parameters, when one considers the caloric curves. Indeed 
one notes in Fig. [7^ that the to = line (full) crosses the magnetized curve (dashed) to the left of T', the temperature 
corresponding to the saddle node bifurcation shown in Fig. ^}p. This leads to the impossibility of applying the usual 
Maxwell construction. However, this is not always the case and, for example, D — 1, r — 1.157 > rtr leads to the usual 
features of a crossing to the right of T' . This confirms that in the interval r g [rtr, 1-157], the mean field cj)^ model has 
a narrow region of negative specific heat, where the transition will be first order in the canonical ensemble and second 
order in the microcanonical. Unfortunately, the points near T' are extremely difficult to obtain because of numerical 
inaccuracies, and we are therefore unable to report a clear determination of negative specific heat. In conclusion, this 
model shows a scenario similar to the BEG model |2lj . in the case of a dynamical Hamiltonian. Meanwhile, similar 
results have been published for some extensions of the HMF model . 



IV. CONCLUSION 



The Blume-Emery-Griffiths mean-field model was shown to be an excellent benchmark to discuss relations between 
canonical and microcanonical ensembles in long range interacting systems '21^ . Indeed, this model is exactly solvable 
in both ensembles and is, at the same time, sufficiently rich to display such interesting features as negative specific heat 
and temperature jumps in the microcanonical ensemble. However, it has no dynamics and only the thermodynamical 




FIG. 7: Thermodynamics of the (j)^ model in the region of the first-order transition (r = 1.18, D — 1.0): Panel (a) presents the 
caloric curve and panel (b) the magnetization as a function of the energy. The critical temperature is Tc = 0.0156 < Ttr- Data 
obtained with microcanonical simulations, A'^ = 1024 t — 1.25 10®. 



behavior can be investigated. This is why we need to study models that displays all these interesting thermodynamical 
effects, but for which one would also dispose of an Hamiltonian dynamics. This point was already addressed in the 
framework of the HMF model and in particular in its two-dimensional version (see ]6| for a review) . Having access to 
dynamics, one can moreover study non equilibrium features. 

The mean-field models that we have considered in this paper are exactly solvable in the canonical ensemble by 
a Hubbard-Stratonovich transformation. The data in the microcanonical ensemble were obtained using molecular 
dynamics simulations. We have shown that these models have first and second order phase transitions and tricritical 
points. Their phase diagram allows to test the presence of ensemble inequivalence near canonical first order phase 
transitions and we have also studied spontaneously generated out-of-equilibrium structures. 
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